##############################
# Figure 1
##############################

rm(list=ls())

#sink("~/Dropbox/COVID Peru/08_replication/04_figures.txt")

library(foreign)
library(Hmisc)
library(ri)
library(stargazer)
library(gridExtra)
library("ggplot2")

##############################
# Prepare data
##############################

# load data
d = read.dta("clean_survey_experiment_covid_peru.dta")

# explore data
names(d)
head(d)
table(d$t_ind_economia)
table(d$t_ind_salud)

############################
# Figures
############################

# plot 1a
acuerdo_cuarentena_binary_e = mean(d$acuerdo_cuarentena_binary[d$t_ind==2], na.rm = T)
acuerdo_cuarentena_binary_e_high = acuerdo_cuarentena_binary_e + 1.96*(sd(d$acuerdo_cuarentena_binary[d$t_ind==2])/sqrt(length(d$acuerdo_cuarentena_binary[d$t_ind==2])))
acuerdo_cuarentena_binary_e_low = acuerdo_cuarentena_binary_e - 1.96*(sd(d$acuerdo_cuarentena_binary[d$t_ind==2])/sqrt(length(d$acuerdo_cuarentena_binary[d$t_ind==2])))
acuerdo_cuarentena_binary_s = mean(d$acuerdo_cuarentena_binary[d$t_ind==1], na.rm = T)
acuerdo_cuarentena_binary_s_high = acuerdo_cuarentena_binary_s + 1.96*(sd(d$acuerdo_cuarentena_binary[d$t_ind==1])/sqrt(length(d$acuerdo_cuarentena_binary[d$t_ind==1])))
acuerdo_cuarentena_binary_s_low = acuerdo_cuarentena_binary_s - 1.96*(sd(d$acuerdo_cuarentena_binary[d$t_ind==1])/sqrt(length(d$acuerdo_cuarentena_binary[d$t_ind==1])))
acuerdo_cuarentena_binary_c = mean(d$acuerdo_cuarentena_binary[d$t_ind==0], na.rm = T)
acuerdo_cuarentena_binary_c_high = acuerdo_cuarentena_binary_c + 1.96*(sd(d$acuerdo_cuarentena_binary[d$t_ind==0])/sqrt(length(d$acuerdo_cuarentena_binary[d$t_ind==0])))
acuerdo_cuarentena_binary_c_low = acuerdo_cuarentena_binary_c - 1.96*(sd(d$acuerdo_cuarentena_binary[d$t_ind==0])/sqrt(length(d$acuerdo_cuarentena_binary[d$t_ind==0])))
acuerdo_cuarentena_binary = c(acuerdo_cuarentena_binary_e,acuerdo_cuarentena_binary_s,acuerdo_cuarentena_binary_c)
acuerdo_cuarentena_binary_high = c(acuerdo_cuarentena_binary_e_high,acuerdo_cuarentena_binary_s_high,acuerdo_cuarentena_binary_c_high)
acuerdo_cuarentena_binary_low = c(acuerdo_cuarentena_binary_e_low,acuerdo_cuarentena_binary_s_low,acuerdo_cuarentena_binary_c_low)
groups = c("Economy","Health","Control")
data = data.frame(acuerdo_cuarentena_binary,acuerdo_cuarentena_binary_high,acuerdo_cuarentena_binary_low,groups)
names(data)
data 

p1a = ggplot(data, aes(groups, acuerdo_cuarentena_binary, fill=as.factor(groups))) +  geom_bar(position="dodge", stat="identity", show.legend = F) + 
                geom_errorbar(aes(ymin=acuerdo_cuarentena_binary_low, ymax=acuerdo_cuarentena_binary_high),
                size=.6,    # Thinner lines
                width=.3,
                position=position_dodge(.9)) + theme_bw() + scale_y_continuous(breaks=c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)) + coord_cartesian(ylim=c(0,1)) + 
theme(legend.key = element_rect(colour = NA))  +  scale_fill_manual(name = "Groups",values=c("grey20", "grey40", "grey60")) +
xlab("") + ylab("Support for quarantine (binary)") + theme(axis.text=element_text(size=18), axis.title=element_text(size=18))
p1a

round(acuerdo_cuarentena_binary_e,3)
round(acuerdo_cuarentena_binary_s,3)
round(acuerdo_cuarentena_binary_c,3)

# plot 1b
acuerdo_cuarentena_e = mean(d$acuerdo_cuarentena[d$t_ind==2], na.rm = T)
acuerdo_cuarentena_e_high = acuerdo_cuarentena_e + 1.96*(sd(d$acuerdo_cuarentena[d$t_ind==2], na.rm = T)/sqrt(length(d$acuerdo_cuarentena[d$t_ind==2])))
acuerdo_cuarentena_e_low = acuerdo_cuarentena_e - 1.96*(sd(d$acuerdo_cuarentena[d$t_ind==2], na.rm = T)/sqrt(length(d$acuerdo_cuarentena[d$t_ind==2])))
acuerdo_cuarentena_s = mean(d$acuerdo_cuarentena[d$t_ind==1], na.rm = T)
acuerdo_cuarentena_s_high = acuerdo_cuarentena_s + 1.96*(sd(d$acuerdo_cuarentena[d$t_ind==1], na.rm = T)/sqrt(length(d$acuerdo_cuarentena[d$t_ind==1])))
acuerdo_cuarentena_s_low = acuerdo_cuarentena_s - 1.96*(sd(d$acuerdo_cuarentena[d$t_ind==1], na.rm = T)/sqrt(length(d$acuerdo_cuarentena[d$t_ind==1])))
acuerdo_cuarentena_c = mean(d$acuerdo_cuarentena[d$t_ind==0], na.rm = T)
acuerdo_cuarentena_c_high = acuerdo_cuarentena_c + 1.96*(sd(d$acuerdo_cuarentena[d$t_ind==0], na.rm = T)/sqrt(length(d$acuerdo_cuarentena[d$t_ind==0])))
acuerdo_cuarentena_c_low = acuerdo_cuarentena_c - 1.96*(sd(d$acuerdo_cuarentena[d$t_ind==0], na.rm = T)/sqrt(length(d$acuerdo_cuarentena[d$t_ind==0])))
acuerdo_cuarentena = c(acuerdo_cuarentena_e,acuerdo_cuarentena_s,acuerdo_cuarentena_c)
acuerdo_cuarentena_high = c(acuerdo_cuarentena_e_high,acuerdo_cuarentena_s_high,acuerdo_cuarentena_c_high)
acuerdo_cuarentena_low = c(acuerdo_cuarentena_e_low,acuerdo_cuarentena_s_low,acuerdo_cuarentena_c_low)
groups = c("Economy","Health","Control")
data = data.frame(acuerdo_cuarentena,acuerdo_cuarentena_high,acuerdo_cuarentena_low,groups)
names(data)
data 

p1b = ggplot(data, aes(groups, acuerdo_cuarentena, fill=as.factor(groups))) +  geom_bar(position="dodge", stat="identity", show.legend = F) + 
  geom_errorbar(aes(ymin=acuerdo_cuarentena_low, ymax=acuerdo_cuarentena_high),
                size=.6,    # Thinner lines
                width=.3,
                position=position_dodge(.9)) +
  theme_bw() + scale_y_continuous(breaks=c(1,2,3,4,5)) + coord_cartesian(ylim=c(0,5)) + 
  theme(legend.key = element_rect(colour = NA))  +  scale_fill_manual(name = "Groups",values=c("grey20", "grey40", "grey60")) +
  xlab("") + ylab("Support for quarantine (ordinal)") + theme(axis.text=element_text(size=18), axis.title=element_text(size=18))
p1b

###############
# Plot 
###############

grid.arrange(p1a,p1b, ncol = 1, top = "")

sink()
